ARC- 147 10 



1 Patent 



ENHANCED ELLIPTIC GRID GENERATION 
Cross Reference To Related Applications : 

This application claims the benefit of U. S. Provisional Application 
60/425,750. 

Origin of the Invention : 

The invention described herein was made by an employee of the United 
States Government and may be manufactured and used by or for the Government 
for governmental purposes without the payment of any royalties thereon or 
therefor. 

Technical Field : 

The present invention is a method for eliminating requirements for 
parameter inputs for generalized grid generation in modeling of engineering 
systems. 

Background of the Invention : 

A large amount of effort has been devoted to developing, enhancing and 
using grid generation techniques, through solution of elliptic partial differential 
equations ("PDEs"). Elliptic grid generation methods generally focus on 
developing body-conforming grids around bodies for simulations of external fluid 
flow. The grids thus generated are smooth, having at least continuous first and 
second derivatives, appropriately stretched or clustered, and are orthogonal over 
most of the grid domain. Inclusion of inhomogeneous terms in the PDEs allows a 
grid to satisfy clustering and orthogonality properties in the vicinity of specific 
surfaces in three dimensions and in the vicinity of specific lines in two dimensions. 

Following the work of Thompson, Thames and Mastin, Jour. Computational 
Physics, vol. 24, 1977, pp. 274-302, three-dimensional governing equations for 
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elliptic grid generation are often expressed as: 

?xx + % y + 5» = Pd,T|,0 = -afsgn(£ - exp{-b i li - y>, (1 A) 

*lxx + *l yy + ^Izz = QCI^.D = -CiSgn(ri - r\) exp{-dilTi - Tljl}, (IB) 
+ £ yy + £ zz = R(I.T|,0 = -cfsgnfc - exp{-f i l£ - y }, (1C) 

where ^, r] and £ are generalized curvilinear coordinates, x, y and z are Cartesian 

coordinates, and P(^,r]£), Q(^,r]X), and R(£,r|£), are inhomogeneous terms, % b i5 

c i5 d b c { and f { are manually selected constants, and the subscript "i" refers to a 

particular boundary component associated with the problem. 

In a two dimensional study by Steger and Sorensen, Jour. Computational 

Physics, vol. 33, 1979, pp 405-410, the authors use the following governing 

equations, 

£x* + % y = -arsgnOi - rii) exp^lr] - r^i}, (ID) 

*lxx + Tlyy = -Ci'Sgn(T| - K}) exp{-d i lTl - 1^1}, (IE) 

for a given r\ boundary. The quantities ^ and Cj are generalized to functions a^) 
and Cj (^), respectively, and the values of these quantities are computed as part of 
the solution by requiring a specified spacing between a given r\ boundary and an 
adjacent grid line, and grid orthogonality at this r\ boundary. In any two 
dimensional problem, the decay parameters, such as d^ must be prescribed or 
manually inserted for each of the boundaries in any coordinate direction. 
However, as will be seen in the subsequent development, these values of d { are 
coupled with the values computed for the quantities % and c { respectively so that 
explicit prescriptions of values for the parameters d { are in conflict with the values 
computed for ^ and for c { . Further, the process of selecting the two values for the 
parameters d { for two opposing boundaries, and, by extension, four values for four 
boundaries in a two dimensional problem, is cumbersome for static grids and is 
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infeasible where dynamically changing grids are required. In three dimensions, the 
parameter values need to be prescribed for six boundaries, one for each of six 
boundaries. 

What is needed is an approach that provides an automatic procedure for 
generating an elliptic grid that does not require manual insertion or user 
prescription of these decay parameters for a two dimensional or three dimensional 
grid generation problem. Preferably, these decay parameters should allow for a 
variable rate of decay from different points on any grid boundary, should arise 
automatically in the formulation and solution of the problem and should permit an 
interpretation in terms of one or more physical quantities associated with the 
problem. 

Summary of the Invention : 

These needs are met by the invention, which provides a process for 
generating an elliptic grid in which (1) grid points tend to cluster near a boundary 
at a desired rate, which may vary from point to point where one or more coordinate 
variables may undergo a relatively large change in value), (2) grid lines 
corresponding to a constant value of a coordinate are approximately parallel to or 
perpendicular to the local boundary line, and (3) user prescription or manual 
insertion of parameters to achieve a desired grid behavior in terms of clustering 
and orthogonality near boundaries is not required (or allowed). The process 
includes the following steps: 

providing defining equations, valid near at least one boundary segment in a 
generalized coordinate system, of a selected grid system, where each of the 
defining equations has at least two independent Cartesian coordinate variables, has 
at least one generalized coordinate as a dependent variable, and comprises a partial 
differential equation, expressed in at least one generalized coordinate; 
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providing a selected group of boundary constraints for the grid system, valid 
near the at least one boundary segment, where a decay parameter for at least one of 
the generalized coordinate dependent variables near the at least one boundary 
segment is determined as part of a solution for the grid system, rather than being 
prescribed initially; 

providing defining equations and selected boundary conditions, having at 
least two independent coordinate variables and at least one dependent variable, for 
steady state heat transfer on a long thin fin, and providing a correspondence 
between the at least two independent coordinate variables for the grid system near 
the at least one grid boundary segment with the at least two independent coordinate 
variables for the heat transfer problem; 

providing a correspondence between a selected power of at least one heat 
transfer coefficient for the heat transfer problem and at least one decay parameter 
for the grid system near the at least one grid boundary segment; and 

determining a solution of the grid system near the at least one grid boundary 
segment that incorporates at least one boundary constraint comprising the at least 
one decay parameter determined for the grid system. 
Brief Description of the Drawings : 

Figure 1 is a flow chart of a procedure for practicing the invention. 

Figures 2 and 3 illustrate grids, computed using the invention, for a two 
dimensional annular region and for a two dimensional convex region, illustrating 
clustering of grid points near an inner boundary and near an upper boundary, 
respectively. 
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Description of Best Modes of the Invention : 

The equations (1A), (IB) and (1C) are modified here in the context of 
equations (ID) and (IE) and are written in the following form, as a six-equation 
set for each of the d,T|,C) boundaries, where a ki = a ki (r|£), ,c ki = c ki d£), and e ki = 
e k i (!,r|) (k = 1, 2, 3). The decay parameters b { , d { and fj are (positive) constants for 
any given boundary segment and are expressed as parameter functions, b^ri,^), 
d,d>£)» and f;d,r|). Without loss of generality, one can assume that near a given |- 
boundary segment i, % - > 0,in a selected region on one side of this boundary 
segment, where sgnd - ^) > 0. Treatment of a situation with sgnd - ^) < 0 is 
analogous. In a region close to this boundary segment, where bj(r|£)li - Q « 1» 
the defining equations and the nonhomogeneous terms have the forms 



ixx + ?yy + ^zz = Pl(^l£), (2A) 

*lxx + llyy + ^lzz = Qld^^X (2B) 

t xx + ^ y + l zz = h(tr\&), (2C) 

p,d,f|.0 = -au(Tl,0 sgnd - li) exp{- bfagfe - 1,1}, (2D) 

a -a 14 (Tj,0 + a„(T|,0 b^) d - li), (2E) 

q,(i,t|.0 = sgnd - li) exp{- Di (Ti£)l£ - (2F) 

* -c u (ti,C) + c u (r,£) b,(Ti,0 d - Si), (2G) 

ri d,ri,U = -e 1(i (Ti,D sgnd -.6) exp{- b^Ol&li - 6 I}, (2H) 

~ -e^) + e 1(i (ri,D b,(T|,0 d - £), (21) 
In a region close to an Y]-boundary segment, where djd,£)lTj - r\) « 1 and r| - t), 
> 0, the defining equations and the nonhomogeneous terms have the forms 

ixx+i yy +izz=p 2 (i>'n£)> (3A) 

Tlxx + "Hyy + Tlzz = q 2 d,Tl,D 5 (3B) 

Cxx + £yy + £zz = tJ&I\®, (3C) 
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P2& T l£) = -a 2li (££) sgn(r] - r\J exp{- d^&m. - r^l}, (3D) 

= + a^fcO d,^,D (n - rfc), (3E) 

fhfon.D = -c 2li (££) sgn(T] - Tii) exp{- di(££)lTi - t^I}, (3F) 

* -c 2ii (|£) + c 2)i (|£) d^Q (n - rii), (3G) 
r 2 (|,Ti,0 = -e 2ii (^) sgn(ri - r^) exp{- d^QIri - t^I}, (3H) 

* -e 2>i (U) + eyd,© di(i,0 (t| - rii), (31) 
In a region close to an ^-boundary segment, where ffar]) \t> - Q « 1 and t, - ^ > 
0, the defining equations and the nonhomogeneous terms have the forms 

ixx + iyy + izz = P3(^'n^), (4A) 

*lxx + "Hyy + *lzz = VsfeW&l (4B) 

£ xx + £ yy + C 2Z = r 3 ^,r|£), (4C) 

P 3 (|,ti,0 = -a 34 (^,Ti) sgnft - Ci) exp{- f,(|,i|)^ - £1}, (4D) 

= -a 3>i (|,ri) + a 3>i (l,r 1 ) f^r]) ft - £), (4E) 

q 3 &ri£) = -c 3ii (i,T,) sgnft - Q exp{- 1^)15 - 5,1}, (4F) 

= -c 3>i (^ri) + c 3>i (^,ri) fi(i,Ti) ft - £), (4G) 

r 3 (i,r)£) = -e 3(i (^,ri) sgnft - Q exp{- f^)^ - y>, (4H) 

« -e 3>i (^n) + e 3>i &ri) f^r]) £ - £), (41) 

Where the preceding approximations, near a boundary segment £ = £j, for 
example, are used, the defining equations for £ become 

£xx + lyy + Izz - H&M) U^v) (£ - Cj) = -a^T)) sgnft - (5A) 

^lxx + *l yy + Tlzz - <Hjg,t\) f&,r}) (£ - 50 = -c 3>i feri) sgnft- Ci), (5B) 

£xx + Cyy + Czz " ffoti ft - = "-^(1*) Sgnft - Q. (5C) 



If one ignores the nonhomogeneous terms on the right hand side of Eq. (5C) 
and ignores the dependence upon z (or £), this relation is seen to represent a steady 
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state heat transfer equation for a long, thin fin of width or height 2L, with 
corresponding "temperature" 9 = ^ - £ l5 discussed, for example, by V. S. Arpaci in 
Conduction Heat Transfer , Addition Wesley, Reading, Mass., 1966, pp. 145-147 
and 201-205: 



where h is a heat transfer coefficient in a selected (z-) direction, k is a thermal 
conductivity coefficient and 5 (« L) is thickness of the fin. Similar equations 
apply for 0 = ^ - and 6 = i] - r\ l9 . 

The heat transfer coefficient h corresponds to or is proportional to, for 
example, a decay parameter, such as the coefficient e 3fi (^,T|)f i (|,ri) in Eq. (41). 
Similar equations are developed for the choices 6 = != - £ r , or 9 = r\ - r\ v 

Where the heat transfer coefficient h in the z-direction is small, the thermal 
gradient is correspondingly large normal to the corresponding boundary in the xy- 
space, which requires a close clustering of isothermal lines, or of the corresponding 
time evolving grid lines. This provides a physical basis for the observation that 
clustering near a given boundary increases as a decay parameter (e.g., m 2 in Eq. 
(7)) decreases, and inversely. However, for a grid whose grid lines show a general 
dependence upon both coordinates, x and y (as in Figures 2 and 3, discussed in the 
following), a similar conclusion is arrived at using the following physical 
argument. As the heat transfer coefficient h in the z-direction decreases, the 
temperature gradient near a y-boundary increases correspondingly, for a given heat 
flux, and this leads to a denser clustering of isothermal lines (or, equivalently, of 
grid lines for the time evolving grid). 

Boundary constraints, valid in a one-sided neighborhood of each of the four 
(two dimensional) or six (three dimensional) boundary segments, are incorporated 



a 2 e/dx 2 + d 2 e% 2 - m 2 e = o, 

m 2 = 2h/kS. 



(6) 



(7) 
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by applying Green's theorem in three dimensions for each of the Eqs. (2E), (3G) 
and (41) for the r| and £ boundaries, respectively: 

J s (dQ/dn) do =/ v {-a u (ri£) sgn(£ - ^) + a 14 (T|,Q b^) 6 } dx, (8) 
/ s (38 /an) da = / v {-c 2>i (|£) sgn(ri - ^ + c 2>i ££) d,(££) 6 } dx, (9) 
J s (de/dn) da = J v {-e 3>i &ri) sgnfc - + c^fon) f^r]) 0}} dx, (10) 
where n refers to a direction that is locally normal to a surface S representing a 
totality of six surfaces including the boundary segments of interest and V is a 
volume enclosed or defined by the totality of these six surfaces. These integral- 
type boundary constraints can be used to calculate the decay parameter analogs, 
a u (T|,£)bi(Ti£), c 2 i (T)£)dj(?£) and e 3 i Cn£)f i (£,ri). When expressed in terms of the 
generalized coordinates (%,r|£), the boundary constraints set forth in Eqs. (8), (9) 
and (10) are transformed, for t, for example, as follows. 

J S I da = J s (d8/dn) da = J s (dtjdn) da (11) 
The integral J s I do in Eq. (14) can be written as an algebraic sum of six 
integrals, evaluated over the indicated boundary segments: 
J S I da = J ?max I da + J nmax I da + J^ I do 

-J^nlda-J^Ida-J^Ida, (12) 
where the surface configurations ^max, §min, etc. represent the corresponding 
boundary segments that together make up the surface S. For the first and fourth 
integral pair, the second and fifth integral pair, and the third and sixth integral pan- 
in Eq. (15), the following respective relations may be verified: 

J^I da = J | (l/jVa n ) a 13 [(„ 2 + y n 2 + z n 2 ) (x^ 2 + y n 2 + z^ 2 )] 1/2 d£ dr] . (13) 
. J, I da = J, (1/Wa 22 ) a 23 [( X| 2 + y 5 2 + z, 2 ) (x K 2 + y ? 2 

+ z^ 2 )] 1/2 d^, (14) 
J^ I da = J t (1/J) [a 33 (x, 2 + y n 2 + z n 2 ) (x^ 2 + y ? 2 + z ? 2 )] 1/2 dr, d| , (15) 
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*n=P<!k 2 + % 2 + V)> (16A) 

a 22 = J 2 (i 1s 2 + T| y 2 + T| I 2 ) > (16B) 

a 33 = J 2 (C x 2 + V + ^ ( 16C ) 

«i2 = J 2 ( + + ^1z)> (16D) 

a 13 = J 2 (U x + ^ y + ? z a (16E) 

<*23 = j2 ( + *l y Cy + (16F) 



where J = J((x,y,z)/$;,r]£)) is a Jacobian of the transformation (x,y,z) — > (^,ip£). 

Equations (13)-(15) can be used to express the boundary constraints in the 
computational space (generalized variables (^,ip£)), analogous to the Eqs. (8)-(10) 
expressed in Cartesian coordinate space. The defining equations in computational 
space that are solved, subject to the boundary constraints set forth in Eq. (41) for a 
^-boundary, are 

a il X i,|§ + a 22 X i,tir] + a 33 X i,K + 2{OC 12 X^ + a 13 X^ + (X 23 X^) = 

= -J 2 { P 3 x i,§ + q 3 Xi.r) + r 3 x iit }, (17) 
Xj = x, yorz. (18) 
Figure 1 is a flow chart of a suitable procedure for practicing the invention 
for a two dimensional or three dimensional grid system. In step 1 1, the system 
provides defining equations, valid near one or more grid boundary segments in a 
generalized coordinate system, of a selected grid system, where each of the 
defining equations has at least two independent Cartesian coordinate variables, has 
at one generalized coordinate as a dependent variable, and comprises a partial 
differential equation, expressed in at least one generalized coordinate. 

In step 13, the system provides a selected group of boundary constraints for 
the grid system, valid near the one or more boundary segments, where a decay 
parameter for at least one of the generalized coordinate dependent variables near 
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the one or more boundary segments is determined as part of a solution of the 
defining equations, rather than being prescribed initially. 

In step 15, the system provides defining equations and selected boundary 
conditions, having at least two independent coordinate variables, for steady state 
heat transfer on a long thin fin, and the system provides a correspondence between 
the at least two independent coordinate variables for the grid system near at least 
one grid boundary segment with the at least two independent coordinate variables 
for the heat transfer problem. 

In step 17, the system provides a correspondence a selected power of at least 
one heat transfer coefficient for the heat transfer problem and at least one decay 
parameter for the grid system near the at least one grid boundary segment. 

In step 19, the system determines a solution of the grid system near the at 
least one grid boundary segment that incorporates at least one boundary constraint 
comprising the at least one decay parameter determined for the grid system. 

Figure 2 illustrates a result of application of the invention to a two 
dimensional annular region to provide a grid in which grid points cluster near an 
inner boundary of the annulus. Where a geometrical system, such as an annulus, 
evolves with time, a grid according to the invention is developed at each of a 
selected sequence of times, with the parameters subject to the boundary constraints 
being allowed to vary from one time to another time. Each of these grids can be 
used to perform a finite element or finite difference analysis on the geometrical 
object that represents the time evolving system at one of these times. 

Figure 3 illustrates a result of application of the invention to a two 
dimensional convex-concave geometry region to provide a grid in which grid 
points cluster near an upper boundary of the region. 
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A decay parameter, such as e 3ii (|,ii)f i (^,Ti), may vary with one or more of 
the generalized coordinates, such as | and/or r\, rather than being constant; and this 
variation is determined as part of the solution of the grid problem, rather than being 
prescribed initially by the user. A grid solution can be determined for a temporally 
constant environment. Alternatively, a time evolving environment can be allowed 
to vary at each of a sequence of times, and a grid solution can be determined for 
each of this sequence of times, using the preceding analysis at each of these times. 

The preceding analysis has focused on neighborhoods of the grid boundary 
segments. As noted in the preceding, in an interior region, far from the grid 
boundary segments, the defining partial differential equations (PDEs) become 
homogeneous, and standard analysis of elliptic PDEs is applied to determine an 
interior solution, which is automatically matched in the solution process across a 
selected interior boundary to the solution obtained for the grid boundary segments. 



